# Load and combine scenario question variables
library(haven)
library(dplyr)
load(file="UMAS0035_IS.RData")

pn$sq <- as.integer(pn$scenario_question_closed)
pn$sq[pn$treatment==6] <- pn$scenario_question_original[pn$treatment==6]

#download replication data for Sagan and Valentino 2017 and merge
#https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/9XHAPW
#This line only needed the first time
download.file("https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/9XHAPW/WOKCHO","Sagan Valentino Replication Data for IS Hiroshima in Iran.dta")
ps <- data.frame(read_dta("Sagan Valentino Replication Data for IS Hiroshima in Iran.dta"))
ps$sq <- 7-ps$prefer_strike_D1_6 #recode because backwards - see S&V notes
ps$weight <- ps$expweight
ps %>% filter(treatment==15) %>% select(sq,weight,caseid,treatment) -> ps
pn <- data.frame(bind_rows(pn,ps)) #merge with prolifnorm replication to compare

#recode
#sq = survey question, ge4 = >=4, i.e., prefer strike
pn$sq.ge4 <- as.integer(pn$sq>=4)

#load analysis libraries
library(survey)

# UMAS0035_IS codebook.docx: variables
# UMAS0035 question ordering.docx: treatments

#Note that all n values reported in the text are raw and do not take into account missing variables.

#H1: Replication (n 55, treatment 15 from Sagan and Valentino, n=160 versus our replication in treatment 6 , n=250)
pn$replicate <- NA
pn$replicate[pn$treatment==6] <- "2018"
pn$replicate[pn$treatment==15] <- "2015"
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(replicate) & !is.na(sq.ge4))
svyttest(sq.ge4~replicate,pn.design,na.rm=TRUE)
svyby(~sq.ge4,~replicate,design=pn.design,na.rm=TRUE,FUN=svymean)

#"Civilians" variant (n 56, treatment 6 versus 1+2)
pn$city.civilians <- NA
pn$city.civilians[pn$treatment==6] <- "City"
pn$city.civilians[pn$treatment<=2] <- "Civilians"
pn.design <- subset(pn.design <- svydesign(ids=~0,data=pn,weights=~weight),!is.na(city.civilians)&!is.na(sq.ge4))
svyby(~sq.ge4,~city.civilians,design=pn.design,na.rm=TRUE,FUN=svymean)
svyttest(sq.ge4~city.civilians,pn.design,na.rm=TRUE)

#H2. Knowledge of legal norm decreases strike support
#80% of those who were asked prior to strike were correct (n 57)
pn$flutest.eq2 <- pn$fluency_test
pn$flutest.eq2[pn$flutest.eq2==1] <- 0 #recode since 2 is the "right" answer (reverse of fluinfo)
pn$flutest.eq2[pn$flutest.eq2==2] <- 1
pn$flutest.eq2.2b2c2d <- pn$flutest.eq2
pn$flutest.eq2.2b2c2d[pn$treatment<8] <- NA #just 8-10 (fluency before sq)
pn.design <- svydesign(ids=~0,data=pn,weights=~weight) #set up survey design
svymean(~pn$flutest.eq2.2b2c2d,design=pn.design,na.rm=TRUE)

#45% of those who answered correctly were in favor (n 57)
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(sq.ge4)&!is.na(flutest.eq2)&pn$treatment<=2)#just 1-2 (fluency after sq)
svyby(~sq.ge4,~flutest.eq2,design=pn.design,na.rm=TRUE,FUN=svymean) #answers to sq conditional on flutest=2
svyttest(sq.ge4~flutest.eq2,design=pn.design)

#H3 - sentiment and survey question - significant difference in means in the latter given the former, i.e., given a "strong" answer, what do they say to the survey?
#iv <- "sentiment"
#dv.nam <- "Agreement that it is wrong to attack civilian populations"
#Sentiment, across treatments 1-9 (n 58)
pn.tab <- as.data.frame(list(sentiment=1:5))
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.tab$svyFreq <- as.vector(svyhist(~sentiment,pn.design,prob=TRUE,breaks=c(0,1,2,3,4,5)+.5)$counts) 
pn.tab$percentage <- pn.tab$svyFreq/sum(pn.tab$svyFreq)
print(pn.tab)

#Sentiment, across treatments 1-6 (n 58)
pn.tab <- as.data.frame(list(sentiment=1:5))
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),treatment<=6)
pn.tab$svyFreq <- as.vector(svyhist(~sentiment,pn.design,prob=TRUE,breaks=c(0,1,2,3,4,5)+.5)$counts) 
pn.tab$percentage <- pn.tab$svyFreq/sum(pn.tab$svyFreq)
print(pn.tab)

#Sentiment, across treatments 7-9 (n 58)
pn.tab <- as.data.frame(list(sentiment=1:5))
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),treatment>=7)
pn.tab$svyFreq <- as.vector(svyhist(~sentiment,pn.design,prob=TRUE,breaks=c(0,1,2,3,4,5)+.5)$counts) 
pn.tab$percentage <- pn.tab$svyFreq/sum(pn.tab$svyFreq)
print(pn.tab)

#H3b - sentiment and survey question - significant difference in means in the latter given the former, i.e., given a "strong" answer, what do they say to the survey?
#39% versus 74% (n 62)
pn$sentiment.eq1 <- pn$sentiment==1
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(sq.ge4)&treatment<=2)
svyby(~sq.ge4,~sentiment.eq1,design=pn.design,na.rm=TRUE,FUN=svymean)
svyttest(sq.ge4~sentiment.eq1,design=pn.design)

#56% in 1-2 agreed w norm (n 62)
pn.tab <- as.data.frame(list(sentiment=1:5))
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),treatment<=2)
pn.tab$svyFreq <- as.vector(svyhist(~sentiment,pn.design,prob=TRUE,breaks=c(0,1,2,3,4,5)+.5)$counts) 
pn.tab$percentage <- pn.tab$svyFreq/sum(pn.tab$svyFreq)
print(pn.tab)

#80 percent of somewhat group support a strike just 1+2 (sentiment after sq) (n 62)
pn$sentiment.eq2 <- pn$sentiment==2
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(sq.ge4)&treatment<=2)
svyby(~sq.ge4,~sentiment.eq2,design=pn.design,na.rm=TRUE,FUN=svymean)

#H4  - n 63
pn$t1a1b.2a <- NA
pn$t1a1b.2a[pn$treatment<=2] <- 1
pn$t1a1b.2a[pn$treatment==7] <- 2
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(sq.ge4)&!is.na(t1a1b.2a))
svyttest(sq.ge4~t1a1b.2a,pn.design,na.rm=TRUE)
svyby(~sq.ge4,~t1a1b.2a,design=pn.design,na.rm=TRUE,FUN=svymean)

pn$t1f.2a <- NA
pn$t1f.2a[pn$treatment==6] <- 1
pn$t1f.2a[pn$treatment==7] <- 2
pn.design <- subset(svydesign(ids=~0,data=pn,weights=~weight),!is.na(sq.ge4)&!is.na(t1f.2a))
svyttest(sq.ge4~t1f.2a,pn.design,na.rm=TRUE)
svyby(~sq.ge4,~t1f.2a,design=pn.design,na.rm=TRUE,FUN=svymean)